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Abstract 

Reduced density matrix functional theory for the case of solids is presented and a new exchange 
correlation functional based on a fractional power of the density matrix is introduced. We show that 
compared to other functionals, this produces more accurate results for finite systems. Moreover, 



it captures the correct band gap behavior for conventional semiconductors as well as strongly 



correlated Mott insulators, where a gap is obtained in absence of any magnetic ordering. 
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One of the most dramatic failures of the usual local density approximation (LDA) or gen- 
eralized gradient type approximations to the exchange-correlation (xc) functional of density 
functional theory (DFT) is the incorrect prediction of a metallic ground state for the strongly 
correlated Mott insulators, of which transition metal oxides (TMOs) may be considered as 
prototypical. For some TMOs (NiO and MnO) spin polarized calculations do show a very 
small band gap (up to 95% smaller than experiments) but only as the result of AFM or- 
dering, however all TMOs are found to be metallic in a spin unpolarized treatment. On 
the other hand, it is well known experimentally that these materials are insulating in na- 
ture even at elevated temperatures (much above the Neel temperature) This indicates 
that the magnetic order is not the driving mechanism for the gap and is just a co-occurring 
phenomenon. A real challenge for any kind of ab-initio theory then is the prediction of an 
insulating state for these strongly correlated materials in the absence of magnetic order. 
Until now the main focus of reduced density matrix functional theory (RDMFT) has been 
on finite systems like atoms and molecules ^, Q, l], Q, 0, 

B, y, h, y, h, 



10 



functionals 

original Miiller functional 
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13. 



14 



11] with various xc 



15] which are essentially modifications of the 



2|. In the present work we extend RDMFT to the case of solid- 
state systems and introduce a new functional which generates not only accurate gaps for 
conventional semiconductors, but demonstrates insulating behavior for Mott-type insulators 
in the non-magnetic phase. 

Formally, the one-body reduced density matrix 7 for a pure state of N electrons is defined 
as (spin degrees of freedom are omitted for simplicity) 



7 (r, r') = N J tt(r, r 2 , . . . , r^V, r 2 , . . . , r N ) d 3 r 2 . . . d 3 r N . (1) 

Diagonalization of this matrix produces a set of natural orbitals 16) , (pi, and occupation 
numbers, n^, leading to the spectral representation 



7 (r,rO=£^<M r M( r ') 



(2) 



where the necessary and sufficient conditions for ensemble iV-representability 17| require 
< rti < 1 for all i, and Y^i n i — N. In terms of 7, the total ground state energy of the 
interacting system is [l^| (atomic units are used throughout) 
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Ml\ = 4- I limVh(r,r')dV+ f p (t)V(t) d 3 r + \ f P ^ )p{ ^ d 3 r d 3 r' + E xc [y], (3) 
Z J r^r' j z J r — r I 

where p(r) = 7(r, r), V is a given external potential and _E XC we call the xc energy functional. 
Minimizing the total energy with -E X c[7] — — § / l7( r ? r ')| 2 /l r ~~ r/ | d 3 r d 3 r' is equivalent to 
the Hartree-Fock (HF) method. The HF functional satisfies the exact condition of the xc 
hole integrating to minus one, however it does not satisfy the condition of convexity, which 

n 

is required by the exact functional |19[|. Miiller proposed |2j a simple alternative to the HF 



functional in which |7(r,r')| 2 is replaced by (7 p (r, r / ))*7 1 p (r,r'), where 7" indicates the 
power used in the operator sense i.e. 



7 (r,r 



(4) 



Miiller's functional satisfies both conditions 



20j for all p with < p < 1. All studies of 



this functional to date have set p = 1/2 (in the rest of the paper we refer to this as Miiller 
functional). As is well known, however, this functional severely over-estimates electron corre- 
lation 0, [a, [l2, 15] and we find that changing the value of p away from 1/2 only exacerbates 
this problem. Therefore in the present work we instead replace |^y(r, r') | 2 by |7 a (r, r')| 2 , 
where 1/2 < a < 1. This simple functional form interpolates between the uncorrelated HF 
limit (a = 1) and the over-correlating Miiller functional (a = 1/2). Remarkably, we find 
that an intermediate value of a exists for which non-magnetic TMO's are insulators as well 
as accurate gaps found for a diverse set of semiconductors and insulators. The price to pay 
for this improved treatment of correlation is that this functional form fails to satisfy the xc 
hole condition exactly 2l|]. However, all numerical tests show that it is at least convex and 
in fact, as shown in Fig. [T]for the case of finite systems (Fig. [T]is a plot for a representative 
finite system, the LiH molecule), the variation of energy as a function of charge obtained 
using current functional (for a = 0.65,0.7) is much closer to the required straight-line be- 



havior 



19j than the Miiller functional or the so called BBC functionals |5j. This type of 



linear behavior has also been recently demonstrated by Cohen et al. 22] for their density 
functional. 

In solids because of the underlying lattice periodicity of the external potential, the one- 
body reduced density matrix for crystals has the symmetry 7(r + T, r' + T) = 7(1-, r'), 
where T is a primitive translation vector, and thus the natural orbitals are also Bloch 
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FIG. 1: (Color online) Energy for LiH molecule as a function of fractional number of excess 

electrons. Results are obtained using various approximations to the xc functional. 

states. Since there are no Kohn- Sham-like equations to solve, a direct minimization over 
these natural orbitals and occupation numbers is required while maintaining the ensemble 
iV-representability conditions. In practical terms, the natural orbitals are expanded in terms 
of a set of previously converged Kohn-Sham states, and optimization of the natural orbitals 
is performed by varying the expansion coefficients. We should stress that the use of Kohn- 
Sham states as a basis is merely a numerical convenience and results do not depend upon 
the starting point in any way. All calculations are performed using the state-of-the-art full- 



potential linearized augmented plane wave (FP-LAPW) method 



the EXCITING code 



24) 



23), implemented within 



For finite systems the chemical potential is defined as [i(N) = dE(N)/dN where E(N) is 
the ground-state energy of N electrons moving in a fixed external potential. For the exact 
functional, fi(N) consists of horizontal lines with possible steps at integer values of N 19]. 
If Nq is the electron number for which the total system is charge neutral, then the step at iVo 
given by E g = lim I/ _ >0 + [M^o + v) ~ M^o ~~ v)} represents the band gap or twice the chemical 
hardness. For periodic solids both E and N are infinite. Moreover, adding a finite charge r\ 
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to each unit-cell while keeping the external potential fixed would lead to an infinitely charged 
system which is not stable. The only accessible quantity is the energy per unit volume, E, 
as a function of excess charge 77 per unit-cell, where E v+ s v (rj) is the self-consistent total 
energy evaluated with an external potential (V) plus the Coulomb potential of a constant 
charge background (5v) which makes the total system charge neutral. Clearly, the quantity 
E and the chemical potential p,(rj) = dEy + $ v {rf) j 'drj are conceptually quite different from the 
corresponding quantities E(rj) and fi(rj) of a finite system because for infinite systems the 
external potential is not kept fixed as 77 is changed. For the case of small excess charge 77, 
the quantity E can be written, to second order in 77, as 

E v+Sv (rj) = E v {<n) + J 5v(v)p N±v (r)d 3 r (5) 

The Coulomb potential from the constant charge (77) background is given by 

Sv(r) =7] J - — ] — d 3 r' (6) 

and Pm(j), the self-consistent M-electron density evaluated with lattice potential plus neu- 
tralizing background potential, can be expressed as 

pN± v {r) = Piv(r) ± ^(r) (7) 

with / n ± (r)d 3 r = 1. Combining Eqs.([5]), ([6]) and ([7]) we arrive at 

E v+ sM = Mv) ± V I ^f^'W . (8) 

J |r — r I 

Derivative of Eq. (jSj) with respect to 77 yields 



. cm for 77 < 

m = p(ri = 0-) + { (9) 

E g + c r rj for 77 > 

with c\ = 2 / n~ (r) / \ r — r'\d 3 rd 3 r' , c r = 2 / n + (r)/|r — r'\d 3 rd 3 r' and /x = dEy/drj. From Eq. 
© we can see that introduction of the background charge has no effect on the size of the 
band gap. Another aspect that is clear from this equation is that unlike for the case of finite 
systems \x does not consist of horizontal lines on either side of 77 = 0. Although Eq. ([9]) is 
only valid for small values of 77, in practice we find a linear behavior for (x{rj) even for larger 
values of 77, a fact which can be attributed to the metallic nature that the system acquires 
on addition or removal of charge. 




FIG. 2: (Color online) Chemical potential (in Ha) versus the excess charge per unit-cell for (a) Cu, 
(b) Si and (2) NiO. The results are obtained using the current functional with a = 0.65. 

In Fig. [2] we show the plot of the chemical potential versus excess electronic charge for 
a prototype metal (Cu), semi-conductor (Si) and Mott insulator (NiO) obtained using the 
current functional with this value of a = 0.65. For Cu one observes a nearly linear behavior 
with a small negative curvature. In contrast, Si and NiO show a qualitatively different 
behavior in which the curvature of the chemical potential changes sign. We interpret this 
as the appearance of a smoothed discontinuity with a linear behavior to the left and to the 
right. This is strikingly different from that of Cu, which remains metallic. This smoothing 
of discontinuities, a consequence of approximating the xc functional, has already been noted 
in work on finite systems [11 ]. Nevertheless, owing to the near linearity of the chemical 
potential on either side of zero excess charge (see Eq. (jHJ)), the fundamental gap may be 
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estimated llj by construction of two tangents as shown in Fig. [21 This allows a rigorous 
test of the functional by comparison with experiment of the fundamental gap of a wide class 
of materials; semiconductors, insulators, and strongly correlated systems. We emphasize 
that for the Miiller and the BBC functionals, without self-interaction correction, all these 
systems turn out to be metallic. 
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FIG. 3: (Color online) Percentage deviation (5) of calculated band gap from experiment. DFT 
results (circles) are with the LDA. RDMFT values are obtained usin g th e current functional with 
a = 0.65,0.7. The GW data (blue vertical lines) is taken from Refs. [2a], a line is drawn between 
the smallest and the largest value. 



We first consider semiconductors and ionic insulators, for which results are shown in 
Fig. EH again using the current functional with a = 0.65. The agreement between the 
calculated gap and experiments is good, with an average deviation of the calculated gap 
9.5%, comparable to the GW method where, if we choose the best (worst) result from the 
literature for each material we find an average deviation of 3% (9%). It must be stressed that 
this good agreement of band gaps is found for a wide range of semiconductors/insulators, 
with gaps ranging from leV (Ge) to 14.2 eV (LiF) and the character of the material changing 
from predominantly covalent to ionic. We find that on increasing the value of a from 0.65 
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to 0.7 the agreement with experiment improves somewhat, with an average deviation of now 
7.5%. In most cases gaps are overestimated with a = 0.7 and underestimated with a = 0.65. 
Reducing the value of a below 0.65 has deleterious consequences; already for a = 0.6 several 
of the semiconductors are metallic. Therefore, there exists a small range of a (0.65-0.7) for 
which the calculated gaps agree well with experiment. 

The TMOs have fundamental gaps ranging from 4.2 eV (NiO) to 2.4 eV (CoO), however 
the origin of the gap is profoundly different from conventional semiconductors/insulators. 



correlations and, towards the end of the 
J]. Gap formation by strong correlations 



Here the gap is opened by strong Mott-Hubbarc 
3c? series, charge transfer effects also contribute 
is largely unaffected by the magnetic state, and so to focus exclusively on Mott-Hubbard 
physics we calculate these materials in the non-magnetic state. As can be seen in Fig. [3] 
the current functional (for both a = 0.65, 0.7) finds TMOs to be insulators. In this case the 
sensitivity of the gap to a is significantly greater, however, as before, the best value of a is 
still in the range 0.65 — 0.7. 

The traditional quasi particle technique used to treat TMOs, the GW method, has only 
been applied to NiO and MnO in presence of AFM ordering. The reason for this lies in 
the fact that all DFT calculations to-date lead to metallic FeO and CoO, and for the GW 
method to produce a gap starting from this metallic ground state would require a diverging 
self-energy. In contrast, RDMFT is a fully ab-initio non-perturbative theory which with 
the current functional and a fixed value of a, seems to capture not only the physics of 
conventional semiconductors/insulators, but also that of strongly correlated materials. What 
still remains is the method for choosing a. It is clear from the present work that the optimal 
value of a lies in a small range (0.65 — 0.70). In order to refine a any further one could do one 
of two things: make a a functional of the density matrix and determine this functional from 
exact properties of £7 X c[7] or fit a to best reproduce results for a wide range of materials. 

It is clear that the current functional improves the value of the band gaps, however it 
is also important that it performs well for properties already adequately described by DFT 
within LDA. An example of such a property is the equilibrium lattice constant. In Table [Dare 
presented the values for this for Si, cubic-BN and diamond. The results for the equilibrium 
lattice parameter obtained with RDMFT, using both Miiller and current functional (a = 
0.7), are as good as the values obtained using DFT. For the current functional with a = 0.65 
the average percentage deviation of calculated results from experiments is 1.2%. 
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Solid 


Expt. 


DFT-LDA Miiller current 


Diamond 


6.74 


6.68 


6.78 6.75 


Si 


10.26 


10.188 


10.49 10.55 


BN 


6.83 


6.758 


6.86 6.82 






1.03 


1.07 1 



TABLE I: Equilibrium lattice parameter (in a.u.). The DFT values are obtained using the LDA. 
RDMFT results are obtained with the Miiller and current functional with a = 0.7. In the last row 
are given the average percentage deviation of calculated results from experiments. 

In conclusion an extension of RDMFT to solids is presented with introduction of a new 
functional. The values of the fundamental gap for various semiconductors and ionic insula- 
tors are dramatically improved from DFT-LDA values and are in as good agreement with 
experiment as GW values. Furthermore, strong Mott-Hubbard correlation is captured in 
this approach, as calculations of the non-magnetic TMOs show. The Miiller functional ap- 
pears to be unsuitable for use with solids, as it fails to produce a gap even for conventional 
semiconductors. 

We have shown that RDMFT is a viable theory for the study of solids where many body 
effects are important, and where DFT based theories have notably failed. This opens up 
many future possibilities, such as the study of High T c superconductors in their under- 
doped Mott insulating phase. It should also, hopefully, stimulate efforts to develop the 
theory formally, including a temperature dependent extension and a method to produce 
quasi particle spectra. 
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